Stability of Multi-hump Optical Solitons 
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We demonstrate that, in contrast with what was previously believed, multi-hump solitary waves 
can be stable. By means of linear stability analysis and numerical simulations, we investigate the sta- 
bility of two- and three-hump solitary waves governed by incoherent beam interaction in a saturable 
medium, providing a theoretical background for the experimental results reported by M. Mitchell, 
M. Segev, and D. Christodoulides [Phys. Rev. Lett. 80, 4657 (1998)]. 
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Self-guided optical beams, or spatial optical solitons, 
are the building blocks of all-optical switching devices 
where light itself guides and steers light without fabri- 
cated waveguides [|l| . In the simplest spatial soli- 
ton is created by one beam of a certain polarization and 
frequency, and it can be viewed as a self-trapped mode of 
an effective waveguide it induces in a medium Q . When 
a spatial soliton is composed of two (or more) modes of 
the induced waveguide ||, its structure becomes rather 
complicated, and the soliton intensity profile may display 
several peaks. Such solitary waves are usually referred to 
as multi-hump solitons; they have been found for various 
nonlinear models of coupled fields ||| . 

In realistic (nonintegrable) physical models, solitary 
waves can become unstable demonstrating self-focusing, 
decay, or a nonlinearity-driven transition to a stable 
state, if the latter exists §. All these scenarios of soliton 
evolution are initiated by exponentially growing pertur- 
bations and they are attributed to linear instability. It 
is usually believed that all types of multi-hump solitary 
waves are linearly unstable, except for the special case of 
neutrally stable solitons in the integrable Manakov model 
|^| . On the contrary, recent experimental results in- 
dicate the possibility of observing stationary structures 
resembling multi-hump solitary waves. This naturally 
poses a question: Were those observations only possible 
because of short propagation distance and a small insta- 
bility growth rate? Definitely, the experimental results 
challenge the conventional view on multi-hump solitary 
waves in different models of nonlinear physics. 

The purpose of this Letter is twofold. First, we study 
the origin of multi-hump solitons supported by incoher- 
ent interaction of two optical beams in a photorefractive 
medium. We find that multi-hump solitons appear via 
bifurcations of one-component solitons and due to the 
process of hump multiplication, when the intensity pro- 
file of a composite soliton changes from single- to multi- 
humped with increasing power. Second, we perform nu- 
merical stability analysis of two- and three-hump solitary 
waves and also find analytically the instability threshold 
for two-hump solitons. We reveal that two-hump soli- 
tary waves are linearly stable in a wide region of their 



existence, whereas all three-hump solitons are linearly 
unstable, and that even linearly stable multi-hump soli- 
tons may not survive collisions. 

In the experiments |?]], spatial multi-hump solitary 
waves were generated by incoherent interaction of two 
optical beams in a biased photorefractive crystal. The 
corresponding model has been derived by Christodoulides 
et al. ||], and it is described by a system of two coupled 
nonlinear equations for the normalized beam envelopes, 
u(x, z) and w(x, z), which for the purpose of our current 
analysis can be written in the following form ||: 
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where the transverse, x, and propagation, z, coordinates 
are measured in the units of (Ld/k) 1 ^ 2 and Ld, respec- 
tively, Ld is a diffraction length, and k is the wavevector 
in the medium. The parameter A is a ratio of the nonlin- 
ear propagation constants, and s is an effective saturation 
parameter. For s — + 0, the system (Q) reduces to the in- 
tegrable Manakov equations Q. 

We look for stationary, z-independent, solutions of 
Eqs. (|l|) with both components u(x) and w(x) real and 
vanishing as |a;| — > oo. Different types of such two- 
component localized solutions, existing for < {A, s} < 
1, can be characterized by the total power, P(X,s) = 
Pi, + Pm, where the partial powers, P u = J_ oa \u\ 2 dx and 
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P w = J_ oo \w\ 2 dx, are integrals of motion. If one of the 
components is small, i.e. w/u ~ e, Eqs. (Q) become de- 
coupled and, in the leading order, the equation for the 
w-component has a solution uq(x) in the form of a fun- 
damental, sec/i-like, soliton with no nodes. The second 
equation can then be considered as an eigenvalue prob- 
lem for the "modes" w n (x) of a waveguide created by 
the soliton uq(x) with the effective refractive index pro- 
file Uq(x)/[1 + suq(x)}. Parameter s determines the total 
number of guided modes and the cut-off value for each 
mode, A„(s). Therefore, a two-component vector soliton 
(uo,w n ) consists of a fundamental soliton and an nth- 
order mode of the waveguide it induces in the medium. 
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Henceforward we denote such a composite solitary wave 
by its "state vector" : |0,n). 
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FIG. 1. Soliton bifurcation diagram for s = 0.8. Horizontal 
line - branch of the fundamental u— soliton. A-B-C - branch 
of |0, 1) solitons. D-E-F - branch of |0, 2) solitons. Inset: 
Transverse profiles of u- (thin), w- (dashed) fields, and total 
intensity (thick), shown for marked points. 

On the -P(A) diagram (for fixed s), continuous branches 
representing |0, rz) solitons emerge at the points of bifur- 
cations A„(s) of one-component solitons (see Fig. 1). It is 
noteworthy that the first-order mode is in fact the lowest 
possible mode of the waveguide induced by the fundamen- 
tal soliton uq(x). This is because the state |0, 0), node- 
less in both components, can exist only in the degenerate 
case A = 1, when Eqs. (]l|) have a family of equal- width 
solutions uq — A(x) sin 9 and wq = A(x) cos 9, with arbi- 
trary 9, and amplitude A satisfying the scalar equation, 
dA/dx = ±s- 1 [log(l + sA 2 ) - s(l - s)A 2 ] 1 / 2 . 

Additionally, indefinitely many families of vector soli- 
tons \m,n), where m =/= n =/= 0, can be formed as bound 
states of phase-locked \0,n) solitons Jlfj| , pl| . Although 
such states do contribute to the rich variety of the multi- 
hump solitons existing in our model, we exclude them 
from our present consideration. 

Families of vector solitons can be found by numerical 
relaxation technique. Some results of our calculations 
are presented in Fig. 1, for |0, 1) and |0, 2} solitons found 
at s — 0.8. Observing the modification of soliton pro- 
files with changing A (see inset in Fig. 1), one can see 
that the modal description of two-component solitons is 
valid only near bifurcation points. For A > A n , the am- 
plitude of an initially small w-componcnt grows and the 
soliton- induced waveguide deforms. It is this purely non- 
linear effect that gives rise to the existence of multi-hump 
solitons. In particular, two- and three-hump solitons are 
members of the soliton families 0, 1) (branch A-B-C) and 
0,2) (branch D-E-F) originating at different bifurcation 
points. At A ~ A n (s), while the w-component remains 
small, all |0, n) solitons are single-humped, as shown in 
Figs. l(a,d). As the amplitude of w grows with increas- 
ing A, the total intensity profile, I(x) = Uq(x) + w^x), 



develops (n + 1) humps [see Figs. l(b,e)], and at suf- 
ficiently large A the u-component itself becomes multi- 
humped [Figs. l(c,f)]. The separation distance between 
the soliton humps tends to infinity as A — > 1. 

To analyze the linear stability of multi-hump solitons, 
we seek solutions of Eqs. (0) in the form of weakly per- 
turbed solitary waves: u(x,z) = Uq{x) + e[F u (x,z) + 
iG u (x, z)] and w(x, z) — w n (x) + e[F w (x, z) 4- iG w (x, z)], 
where e <C 1. Setting F u ^ w ~ f u ,w{x)e l3z , G u>w ~ 
9u,w{x)e P z , one can obtain the following eigenvalue prob- 
lem (EVP) 
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where a = c = 1/ (1 + si), b = 0, a± = a + 2«„/ (1 + 
si) 2 , ci = c + 2w 2 J{l+sI) 2 , andbi = -2u Q w n /(l + sI) 2 . 
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FIG. 2. Existence and stability domains for two- and 
three- hump solitons. Shown are the existence thresholds Ai (s) 
and A2(s) for the |0, 1) and |0, 2) soliton families. Dashed - 
the line where |0, 1) solitons become two-humped. Shaded 
- analytically obtained instability domain for two-hump soli- 
tons. Squares and circles- numerically obtained instability 
thresholds for |0, 1) and |0, 2} solitons, respectively. 

Because CiCq and LqCi are adjoint operators with 
identical spectra, we can consider the spectrum of only 
one of these operators, e.g. CiCq. Considering the com- 
plex A-plane, it is straightforward to show that A 6 
(-co, —A 2 ) is a continuum part of the spectrum with 
unbounded cigenfunctions. Stable bounded eigenmodes 
of the discrete spectrum (the so-called soliton internal 
modes |l2j|) can have eigenvalues only inside the gap, 
—A 2 < A < 0. The presence of either positive or complex 
A implies soliton instability, because in this case there 
always exists at least one eigenvalue of the soliton spec- 
trum with Re/? > 0. 
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Numerical solution of the EVP (g) shows that both 
1 0, 1 ) and 1 0, 2) types of solitary wave solutions can be 
stable in a certain region of their existence domain, see 
Fig. 2. In the case of |0, 1} solitons, the appearance 
of the instability is related to the fact that close to the 
curve where the total intensity / becomes two-humped 
[dashed line in Fig. 2], a pair of internal modes split 
from the continuum into the gap. As A grows, the corre- 
sponding, purely imaginary, eigenvalues (3 = ±iy/\A(X)\ 
tend to zero, and at a certain critical value A = A cr (s), 
they coincide at f3 = 0. At this point, an eigenmode 
with positive eigenvalue A emerges, thus generating lin- 
ear instability (see Fig. 3) with the instability growth 
rate /3 = a/A(A). For |0, 2) solutions, the dynamics of 
internal modes can not be related in any obvious way 
with a change in the spatial solitary profiles, neverthe- 
less the scenario of the instability development is similar 
to that for two-hump solitons. The dependence of (3 on 
A, for |0, 1} and |0, 2) soliton families giving rise to two- 
and three-hump solitary waves, is shown in Fig. 3 for 
s = 0.3 and s — 0.8, respectively. A decline in the insta- 
bility growth rate as A — > 1 (see Fig. 3) is caused by the 
fact that, in this limit, all multi-hump solitons decompose 
into a number of the neutrally stable |0, 0) solitons sepa- 
rated by infinitely growing distance. Numerical analysis 
in the close vicinity of this limit is unfeasible due to lack 
of computational accuracy. 
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FIG. 3. Instability eigenvalues vs. 
solitons; dashed line - Im/3, bold line ■ 



A for |0, 1) and |0, 2} 
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Note that, within the gap of the continuous spectrum, 
there exist several soliton internal modes not participat- 
ing in the development of the linear instability. Analysis 
of their origin and influence on the soliton dynamics is 
beyond the scope of the present Letter. 

With the aid of analytical asymptotic technique Q , it 
is possible to show that a perturbation mode with small 
but positive eigenvalue, and therefore the linear instabil- 
ity of a general localized solution (u, w), appears if the 
functional J(u,w), defined as 
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changes its sign. The threshold condition J = is, in 



fact, the Vakhitov-Kolokolov stability criterion pj), gen- 
eralized for the case of two-parameter vector solitons. In 
this case, it does not necessarily give a threshold of lead- 
ing instability |T^ |. Therefore, the presence of other in- 
stabilities (which are not associated with the condition 
J = and can have stronger growth rates) is still possi- 



ble, as in some other cases |16 



For two-hump solitons, we have been able to locate the 
critical curve in (A, s)-plane corresponding to the condi- 
tion J = 0. Superimposing this curve onto the numeri- 
cally calculated values A cr (s), we have found a remarkable 
agreement between the numerical and analytical instabil- 
ity thresholds, as shown in Fig. 2. This gives us the first 
example of the generalized Vakhitov-Kolokolov criterion 
for the instability threshold of vector multi-hump solitary 
waves. For the whole family of |0, 2) solutions, including 
three-hump solitons, it appears that J =/= throughout 
the entire existence region. Thus, appearance of insta- 
bility of three-hump solutions is not associated with the 
change of the sign of the functional J. 



To analyze long-term evolution of multi-hump solitary 
waves, we perform numerical simulations of the beam 
propagation for |0, n) solitons within the existence do- 
main A n < A < 1, at fixed s. First, we use no pertur- 
bation so that the soliton instability can only develop 
from numerical noise. As long as the soliton maintains 
its single-humped shape [see corresponding profiles in 
Fig. l(a,d)], it remains almost insensitive to numerical 
noise. Moreover, while the |0, 1} solitons do become two- 
humped at A < A cr , they still remain stable in a wide 
domain of their parameters until the linear instability 
threshold is reached. On the contrary, |0, 2) solitons re- 
main single-humped up to the instability threshold value 
A = A cr , so that all three-hump solitons are indeed unsta- 
ble. Above the instability threshold (i.e. for A cr < A < 1), 
a two-hump soliton splits into two independent single- 
humped beams as a result of the instability developed 
from noise [see Fig. 4(a)], whereas a three- hump soliton 
exhibits a more complex symmetry-breaking instability, 
as shown in Fig. 4(b). 




FIG. 4. Noise-induced splitting of (a) two-hump soliton 
[Fig. 1(c)] and (b) three-hump soliton [Fig. 1(f)]. 
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FIG. 5. Collisions of (a) linearly stable 0, 1) solitons at 
s — 0.8, A = 0.72 with the relative transverse velocity 
v — 0.05, and (b) linearly stable |0, 2) solitons at s = 0.8, 
A = 0.35, with v = 0.09. 

Next, we propagate two-hump (at s = 0.3) and three- 
hump (at s = 0.8) solitons perturbed by an eigenmode 
with the largest instability growth rates, i.e. /3 ma x ~ 
0.055 and /3 max ~ 0.153, respectively. We find that 
in the presence of ~ 6% amplitude perturbation, the 
diffraction-induced decay of a soliton can be stabilized 
by the nonlinearity, whereas its splitting is significantly 
speeded up by the perturbation, compared with splitting 
due to a numerical noise. 

To make a link between our stability analysis and ex- 
periment, we note that for the experiment the diffrac- 
tion length is defined as Ld — 2/sb and nonlinearity of 
the medium (SBN:60 crystal) is characterised by the pa- 
rameter b — kr c ffn^Eo, where r e s is the effective electro- 
optic coefficient (= 280 pm/V), n\> is the background 
refractive index (= 2.3), and E is the applied electric 
field (« 2xl0 5 V/m). For strong saturation we have 
s ~ 1 and Ld « 0.2mm. Now, the characteristic in- 
stability length z cr can be defined through the maximum 
growth rate /3 ma x and, as a result, for two- hump solitons 
at s = 0.3 we obtain z cr « 12.18 mm. These estimates 
indicate that the instability, if it exists, could be detected 
for two-hump solitons within the experimental setup of 
Ref. Q and therefore stable two-hump solitons have been 
indeed observed. 

Importantly, three-hump solitons so far generated in 
the experiment belong to a different class of vector soli- 
tons which, in our notation, can be identified as |1,2) 
states. The extensive numerical analysis of soliton states 
|1, 2) |llj] shows that all such solitons are linearly unsta- 
ble. However, the observation of this instability is beyond 
the experimental parameters of Ref. @] . 

The complex structure of multi-hump solitons and non- 
integrability of the model ([!]) result in a variety of colli- 
sion scenarios, which are quite dissimilar to the collisions 
of multi-hump solitons of the exactly integrable Manakov 
system || . For instance, even linearly stable vector soli- 
tons do not necessarily survive soliton collisions. In Figs. 
5(a,b) we show two examples of non-elastic interaction of 
linearly stable |0, 1) and 0,2) solitons. 

In conclusion, we have analyzed, analytically and nu- 
merically, the stability of multi-hump optical solitons in a 
saturable nonlinear medium. We have found that multi- 



hump solitons are members of an extended class of vector 
solitons which can be linearly stable in a wide region of 
their existence, although they may be destroyed in col- 
lisions. We believe that this is an important physical 
result that calls for a revision of our understanding of 
the structure and stability of many types of multi-hump 
solitary waves in nonintegrable multi-component models, 
usually omitted in the analysis because of their a priori 
assumed instability. 
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